In [1]:
import numpy as np
import plotly.plotly as py
import plotly.offline as po
import plotly.graph_objs as go
po.init_notebook_mode(connected=True)
In an exploration of the possibilities of Jupyter Notebooks, I like to consider NURBS. These curves (or surfaces or solids) play a basic but important role for my graduation thesis. My goal is to make extensive use of numpy and its vectorization capabilities.
NURBS are formed as a combination of basis functions, each defined by the Cox-De Boor recursion formula[1][2].
$$ N_{i,n}\left( u\right) =\frac{u-k_i}{k_{i+n}-k_i}N_{i,n-1}\left( u\right) +\frac{k_{i+n+1}- u}{k_{i+n+1}-k_{i+1}}N_{i+1,n-1}\left( u\right) $$This formula yields the $i$th basis function of polynomial order $n$, given a knot vector $k$. An implementation for arrays of $u$ values is suitable for plotting purposes.
In [2]:
# Cox-De Boor recursion formula, result evaluated for array u
# i = 0, 1, 2, ..., n
def cox_de_boor(knots, i, n, u):
if n == 0:
isbetween = np.logical_and(u >= knots[i], u < knots[i + 1])
return np.where(isbetween, 1, 0)
result = np.zeros_like(u)
# determine if denominators are zero first
Ad = knots[i + n] - knots[i]
if Ad != 0:
A = u - knots[i] / Ad
result += A * cox_de_boor(knots, i, n - 1, u)
Bd = knots[i + n + 1] - knots[i + 1]
if Bd != 0:
B = knots[i + n + 1] - u / Bd
result += B * cox_de_boor(knots, i + 1, n - 1, u)
return result
Constructing NURBS requires taking a sum of the basis functions belonging to a certain polynomial order. It is therefore often convenient to collect values of all basis functions for all $u$ positions in a 2D matrix.
In [3]:
def basis_matrix(knots, u):
# determine polynomial order from assumed open knot vector
n = np.sum(knots == knots[0]) - 1
# number of basis functions
nbas = knots.size - n - 1
result = np.empty((nbas, u.size))
for i in range(nbas):
result[i, :] = cox_de_boor(knots, i, n, u)
return result
Plotting a NURBS basis then becomes straightforward.
In [4]:
knots = np.array([0., 0., 0., 0., 1., 1., 1., 1.])
u = np.linspace(knots[0], 0.99999 * knots[-1], 100)
basis = basis_matrix(knots, u)
data = []
for i in range(4):
data.append(go.Scatter(
#x = random_x,
y = basis[i,:],
mode = 'lines',
name = f'basis function {i}'
))
po.iplot(data, filename='basis')
The actual NURBS is formed as a linear combination of the basis. In a way, each basis value is 'weighted' and combined with a control point.[1]
$$ C\left(u\right) =\sum\limits^{n+1}_{i=1}\frac{N_{i,n}\left(u\right) w_i}{\sum^{n+1}_{j=1}N_{j,n}\left(u\right) w_j}\boldsymbol{P}_i $$
In [5]:
class NURBS(object):
def __init__(self, knots, cpoints, weights):
self.knots = knots
self.cpoints = cpoints
self.weights = weights
def eval(self, u):
basis = basis_matrix(self.knots, u)
wbasis = np.einsum('i, ik -> ik', self.weights, basis)
W = wbasis.sum(axis=0)
return np.einsum('ij, ik -> jk', self.cpoints, wbasis / W)
Note the use of the einsum() function.[3] While this method may require some work to understand, the author believes it provides the computation with more transparency than conventional matrix manipulation. Compared to the formula however, one needs extra indices because of handling many $u$ values at once.
This NURBS object may then be used to plot a curve in space. A spiral for instance, is readily constructed.
In [6]:
knots = np.array([0., 0., 0., 1., 1., 2., 2., 3., 3., 4., 4., 5., 5., 6., 6., 6.])
cpoints = np.array([[3, 0, 0], [3, 3, 1], [0, 3, 2], [-3, 3, 3],
[-3, 0, 4], [-3, -3, 5], [0, -3, 6], [3, -3, 7],
[3, 0, 8], [3, 3, 9], [0, 3, 10], [-3, 3, 11],
[-3, 0, 12]])
w = 1 / np.sqrt(2)
weights = np.array([1, w, 1, w, 1, w, 1, w, 1, w, 1, w, 1])
curve = NURBS(knots, cpoints, weights)
u = np.linspace(knots[0], 0.99999 * knots[-1], 100)
x, y, z = curve.eval(u)
trace = go.Scatter3d(
x = x,
y = y,
z = z,
mode = 'lines',
name = 'NURBS'
)
data = [trace]
po.iplot(data, filename='nurbs')
NURBS may be expanded to higher dimensional objects[1]. The basis is then a 'tensor product' of one dimensional basis functions. For this reason, NURBS surfaces (and solids) are said to have a 'tensor product structure'. As an example, a NURBS surface is computed.
$$ S\left( u,v\right) =\sum\limits^{n+1}_{i=1}\sum\limits^{m+1}_{j=1}\frac{N_{i,n}\left( u\right) N_{j,m}\left( v\right) w_{i,j}}{\sum^{n+1}_{p=1}\sum^{m+1}_{q=1}N_{p,n}\left( u\right) N_{q,m}\left( v\right) w_{p,q}}\boldsymbol{P}_{i,j} $$
In [7]:
class Surface(object):
def __init__(self, knots_u, knots_v, cpoints, weights):
self.knots_u = knots_u
self.knots_v = knots_v
self.cpoints = cpoints
self.weights = weights
def eval(self, u, v):
basis_u = basis_matrix(knots_u, u)
basis_v = basis_matrix(knots_v, v)
wbasis = np.einsum('ij, ik, jl -> ijkl', self.weights, basis_u, basis_v)
W = wbasis.sum(axis=(0, 1))
return np.einsum('ijk, ijlm -> klm', self.cpoints, wbasis / W)
In [8]:
knots_u = np.array([0., 0., 0., 1., 1., 2., 2., 2.])
knots_v = knots_u
cpoints = np.array([
[[0, 0, 0] , [0, 0, -1], [1, 0, -1], [2, 0, -1], [2, 0, 0]],
[[0, 1, -0.5], [0, 1, -1], [1, 1, -1], [2, 1, -1], [2, 1, -0.5]],
[[0, 2, 0] , [0, 2, 0] , [1, 2, 0] , [2, 2, 0] , [2, 2, 0]],
[[0, 3, 0.5] , [0, 3, 1] , [1, 3, 1] , [2, 3, 1] , [2, 3, 0.5]],
[[0, 4, 0] , [0, 4, 1] , [1, 4, 1] , [2, 4, 1] , [2, 4, 0]],
])
w = 1 / np.sqrt(2)
weights = np.array([
[1, w, 1, w, 1],
[1, w, 1, w, 1],
[1, w, 1, w, 1],
[1, w, 1, w, 1],
[1, w, 1, w, 1],
])
surface = Surface(knots_u, knots_v, cpoints, weights)
u = np.linspace(knots_u[0], 0.99999 * knots_u[-1], 100)
v = np.linspace(knots_v[0], 0.99999 * knots_v[-1], 100)
x, y, z = surface.eval(u, v)
trace = go.Surface(
x = x,
y = y,
z = z,
name = 'NURBS surface'
)
data = [trace]
po.iplot(data, filename='surface')